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In 2002 Biskup et al. [Europhys. Lett. 60, 21 (2002)] sketched a rigorous proof for the behavior 
of the 2D Ising lattice gas, which is equivalent to the ordinary spin-1/2 Ising model, at a finite 
volume and a fixed excess SM of particles (spins) above the ambient gas density (spontaneous 
magnetisation). By identifying a dimensionless parameter A(8M) and a universal constant A c , they 
showed in the limit of large system sizes that for A < A c the excess is absorbed in the background 
("evaporated" system), while for A > A c a droplet of the dense phase occurs ("condensed" system). 
By minimising the free energy of the system they derive an explicit formula for the fraction A(A) 
of excess particles forming the droplet. 

To check the applicability of the analytical results to much smaller, practically accessible system 
sizes, we performed several Monte Carlo simulations for the 2D Ising model with nearest-neighbour 
couplings on a square lattice at fixed magnetisation M. Thereby, we measured the largest minority 
droplet, corresponding to the condensed phase, at various system sizes (L = 40, . . . , 640). With an- 
alytic values for for the spontaneous magnetisation mo, the susceptibility \ and the Wulff interfacial 
free energy density tw for the infinite system, we were able to determine A numerically in very good 
agreement with the theoretical prediction. 

Furthermore, we did simulations for the spin-1/2 Ising model on a triangular lattice and with 
next-nearest-neighbour couplings on a square lattice. Again, finding a very good agreement with the 
analytic formula, we demonstrate the universal aspects of the theory with respect to the underlying 
lattice. For the case of the next-nearest-neighbour model, where rw is unknown analytically, we 
present different methods to obtain it numerically by fitting to the distribution of the magnetisation 
density P(m). 



PACS numbers: 05.70.Fh,02.70.Uu,75.10.Hk 

I. INTRODUCTION 

The formation and dissolution of equilibrium droplets 
at a first-order phase transition is one of the longstand- 
ing problems in statistical mechanics Quantities of 
particular interest are the size and free energy of a "crit- 
ical droplet" that needs to be formed before the decay 
of the metastable state via homogeneous nucleation can 
start. For large but finite systems, this is signalised by 
a cusp in the probability density of the order parame- 
ter <j) towards the phase-coexistence region as depicted 
in Figs. Q] and [5] for the example of the two-dimensional 
(2D) Ising model, where <j) = m is the magnetisation. 
This evaporation/condensation "transition point" sepa- 
rates an "evaporated" phase with many very small bub- 
bles of the "wrong" phase around the peak at (£>o from 
the "condensed phase" phase, in which a large droplet 
has formed; for configuration snapshots see Fig. [31 The 
droplet eventually grows further towards = until 
it percolates the finite system in another droplet /strip 
"transition" . The latter transition is indicated in the 2D 
Ising model by the cusp at the beginning of the flat two- 
phase region around m = (see Fig. [T|) . 

Building on the seminal work by Fisher [l| developing 
the droplet picture, early numerical studies of the evapo- 
ration/condensation transition by Binder, Kalos and Fu- 
rukawa date back to the beginning of the 1980s. Re- 
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FIG. 1: Schematic plot of the probability density P(m) of the 
magnetisation in logarithmic form. The marked box indicates 
the position of the cut-out displayed in Fig. [2] The vertical 
(green) line indicates the droplet/strip transition point for 
positive magnetisation m > 0, the use of which will be ex- 
plained later on in Sec. IIII CI 



cently this problem has been taken up again by Neuhaus 
and Hager [4J who discussed it with emphasis on possi- 
ble Gibbs-Thomson and Tolman corrections. This stimu- 
lated further new theoretical [ESQ and numerical [1, 
work. 
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FIG. 2: Probability density of the magnetisation for the two- 
dimensional Ising model around its right peak for different 
system sizes L at the temperature T = 1 .5. The cusp indi- 
cates the evaporation/condensation transition region. On the 
right side of the cusp (evaporated system) a Gaussian peak 
is clearly visible, while on the left side (condensed system) a 
stretched exponential behavior can be seen. The two arrows 
on the x-axis indicate for L = 640 the range of data points 
shown in Fig. 1141 



Here, we follow the exposition of Biskup et al. [H, Q , 
who present their results both in a phenomenological 
liquid-vapour (or solid-gas) picture and also explicitly in 
terms of the simple Ising (lattice-gas) model. The dis- 
tinguishing feature of their work is the formulation of 
a proper equilibrium theory which does not need to ex- 
plicitly involve correction effects a la Gibbs-Thomson or 
Tolman as was done in earlier works [ll|, [H, [l3| • We 
consider this feature as one of the main merits of their 
formulation which can be shown to be equivalent (at least 
in leading order) to the earlier less rigorous treatment in 
1- 

The price one has to pay, however, is a rather intri- 
cate rescaling of the original problem which requires in 
numerical work great care with details. To set the theo- 
retical grounds for our Monte Carlo simulation study and 
in particular to develop intuition for the final representa- 
tion of our results in Figs. IT^HTTl we therefore start first 
with a brief summary of the Biskup et al. [H, 0] theory. 
In order to do so, we restrict ourselves to the special case 
of the 2D Ising model with Hamiltonian 



H = -jJ2s i s j , (1) 



where = ±1 and denotes a (next-)nearest- 

neighbour pair. If a down-spin (<Xj = —1) is treated as a 
particle and an up-spin (<7j = 1) as a vacancy, the system 
can be interpreted as a lattice gas of atoms. 



FIG. 3: Two snapshots of a 320 x 320 n.n. Ising system at 
T = 1.5 and the same value of the magnetisation m = 0.9801 
chosen to be in the vicinity of the evaporation/condensation 
point. Left: Evaporated system, a large number of very small 
excitations (bubbles) exist (1 to 4 spins) and the largest clus- 
ter consists of 5 connected spins. Right: Condensed system, a 
single large droplet with volume 400 spins that has absorbed 
a large amount of the small bubbles. 

II. THEORY 

In this section we summarise the considerations of 
Biskup et al. [j| but specialised for the case of the two- 
dimensional Ising model (not necessarily on a square lat- 
tice). 

We image the following situation: an unconstrained 
[43| Ising system of size V = L x L in the low-temperature 
phase at the inverse temperature (3 = J/k^T > j3 c . If the 
majority of spins is positive (<jj = 1), i.e., the system is 
in the phase with positive magnetisation, then, due to 
thermal fluctuations, there are always some overturned 
negative spins and the total magnetisation is M = itiqV, 
with mo < 1. Here, too = roo(/3) > denotes the 
infinite-volume equilibrium magnetisation (spontaneous 
magnetisation) as, e.g., calculated analytically by On- 
sager and Yang for the square lattice with next-neighbour 
interactions (see Sec. Now, if some volume vl of the 
systems is inverted [44j |. then the magnetisation of this 
constrained system is 

M = m {V - v L ) ~ m v L ■ (2) 

It is important to note, that here we did not require the 
inverted volume vl to be connected or to be of the form 
of a droplet. Still, we present in Fig. 2] this extreme case 
to make it simpler to identify the quantities introduced 
here. Secondly, as only spins are inverted, by symme- 
try it must hold exactly to ' = — 771q , otherwise, an 
completely inverted system would have another value for 
the spontaneous magnetisation. Now, the difference to 
the original, unconstrained system with magnetisation 
Mq = moV is 

6M = M — M a = -2v L m . (3) 

The factor 2 is due to the definition of the Ising spins, 
having a value ±1. The interpretation of this formula is 
as follows: a system which has a difference in the mag- 
netisation of SM to an unconstrained Ising system has a 
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FIG. 4: Ising system of size V with a minority droplet of 
volume dl of negative spins surrounded by positive spins with 
a volume (V — vl), shown in the extreme case where the total 
excess in magnetisation is concentrated in the droplet, i.e. 

Vd = V L ■ 



volume vl of inverted spins. Biskup et al. show that for a 
given magnetisation M the total volume of inverted spins 
Vl can be (in the thermodynamic limit of large systems) 
divided into two parts, unconnected small fluctuations 
with volume Vf and a single large connected droplet with 
volume fd, since there exist no droplets of intermediate 
size @. For the total volume of inverted spins holds 

V L = Vf + v d - 

Now, the free energy can be decomposed according to 
the two contributions. For the droplet it is written as 



Fa 



(4) 



where 7\y is the interfacial free energy per unit volume of 
an ideally shaped droplet, also known as the free energy 
of a droplet of Wulff shape [3] . The contribution of the 
fluctuations is derived in the following manner. From the 
volume V of the whole system already Vd is occupied by 
the single large droplet. The rest of the system has an 
unconstrained magnetisation of Mq = (V — tid)roo< If 
some volume Vf of the remaining spins is inverted, then 
the magnetisation is 



M = (V — Vd — Vi)mo — moVf 



(5) 



Then, the difference SM l to the unconstrained magneti- 



sation Mq is 



(0) 



The contribution to the free energy due to these fluctua- 
tions can be written as 



F { = 



(M l -Ml) 2 _ 2m 2 vf 



2 X V 



(7) 



where x = x(fi) ~ PV[(m?) — (m) 2 ] is the susceptibility 
in the thermodynamic limit. 

Now, the relative volume of the droplet compared to 
the total volume of overturned spins vl is defined as 



A = — or 

VL 



v d = \v L 



(8) 



Hence, Vf can be written as 



v f = v L - v d 



VL 



I'd 



1 - — = v L (1 

VL, 



Using this relation, the total free energy F 



F =T W y/Vd~ + 



2m\v\ 
XV 



=rwVA^+^!(l-A) 2 

or, in the form of Biskup et al., 

F(X) =t w ^I0a(A) 



with 



and 



A (A) = \/A + A (1 - A) 2 



A) . (9) 
F d + F f is 
(10) 

(11) 

(12) 
(13) 



A = 



o™2 3/2 
2 WqV 

xVt w ^/vJ xVtw 



2mlv\ 



(14) 



Now, if the magnetisation is fixed to some value, then the 
total number of overturned spins is also fixed and using 
Eq. © it holds 



VL 



V 



M 
to 



(15) 



As mo, x an d Tw are constants, the only varying quantity 
in Eq. (fT2")) is the relative volume of the droplet A. A 
fully equilibrated thermodynamic system always stays in 
the minimum of the free energy. Therefore, the physical 
Aa, i.e., the correct distribution of overturned volume 
between the droplet and the fluctuations, minimises F in 
the range A € [0,1]. Consequently, the solution of this 
problem is either given by ^i^- — 0, which is 



1 



2\/A 



2A(1 - A) = 



(16) 



or it is one of the boundary values 0, 1. Solving Eq. |T| 
shows that for A < A c the correct solution is A = 0, 
i.e., pure fluctuations and no droplet at all. The point 
A c it given by the condition </>a c (0) = <^A C (A C ) which is 
A c = \/A7+A c (l-A c ) 2 or 



1 



OT c (2 - Ac) 



(17) 



This can be substituted in Eq. (|16[) resulting in 

2(1-A C ) 

VaT(2-a c ) 



2v / a7 



or 



A r 



Inserting this value into Eq. (|T7|) gives 
A c = W^ = 0.918558.. 



(18) 



(19) 
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FIG. 5: Fraction of the excess magnetisation in the largest 
droplet A in dependence of the scaling parameter A. For A < 
A c there is no largest droplet, only fluctuations. At A = A c 
a droplet is formed, containing 2/3 of the total excess. In the 
case A > A c the fraction of the excess is given by Eq. (|20|l . 
The (blue) lines approaching A for A > A c are the Taylor 
series of Eq. (|20|l up to order 4 around A = oo that have the 
form A = 1 - 1/4A - 1/32A 2 - 5/512A 3 - 1/256A 4 + . . . . 



For A > A c the solution is 

X 4 2 

A = - cos 



vr-cos-if^) 



(20) 



These results give rise to the following physical picture. 
For fixed, magnetisation M w Mq, where A(M) < A c , 
the systems contains no droplet, only fluctuations are 
present. At some value M c with A(M C ) = A c two states 
coexist, the state of pure fluctuations and a mixed state 
composed of a droplet that absorbs 2/3 of the fluctu- 
ations and the remaining 1/3 of the fluctuations. For 
smaller magnetisation, i.e. A(M) > A c , the droplet 
grows and thereby absorbs more and more of the back- 
ground fluctuations. The predicted behavior of A = A(A) 
is shown in Fig. [5] 



III. SET UP 

In this work we wanted to answer two questions. On 
the one hand, we wanted to test from which system sizes 
on the theoretical results presented in the last section 
start to yield a good description of the data for the two- 
dimensional Ising model. On the other hand, we wanted 
to check the universal aspects of the theory by using 
different lattice models, namely the triangular nearest- 
neighbour (n.n.) lattice and the next-nearest neighbour 
(n.n.n.) square lattice. In order to do so, A, the fraction 
of the excess of magnetisation in the largest droplet de- 
fined in Eq. (JSj> , had to be measured in dependence of 
the parameter A defined in Eq. (fT4J) . 

To get the correct scaling for the abscissa, the param- 
eter A(vli m 0i X) T w) had to be calculated according to 



Eq. ([14]) . While vl is a free parameter, the magnetisa- 
tion, the susceptibility and the free energy of the Wulff 
droplet per unit volume must be obtained analytically or 
by other means, e.g., as results of simulations. For the 
free ener gy of the Wulff droplet the analytic expression 
E w = 2\/W"S, e.g. [H, can be used. Here, E is the 
volume of the droplet and W is the volume bounded by 
the Wulff plot. Putting E = 1 gives the interfacial free 
energy per unit volume 

7w(/9) = 2VW. (21) 

In the following three subsections we discuss for the 
three studied models the origin of the constants in 
question. For the standard Ising model with nearest- 
neighbour couplings on a square lattice and the Ising 
model on a triangular lattice all relevant constants are 
known from literature, either analytically or from quite 
long series expansions. This is not the case, however, 
for the n.n.n. Ising model and, therefore, here we had to 
apply simulations to retrieve the values. 



A. Parameters for the n.n. Ising model on a square 

lattice 



The critical temperature of the Ising model was given 
in 1941 by Kramers and Wannier [17fl . Using self-duality 
arguments they obtained the expression 

(22) 



hl(l + y/2) 



For the spontaneous magnetisation mo there exists the 
famous Onsager-Yang analytic solution [H, [l^] 



mo(/3) = [1 - suffr 4 (2/?)] 



1/8 



(23) 



Also the susceptibility \ is virtually known to arbitrary 
precision from very long series expansions, e.g., Orrick et 
al. [2(| give the formula 



X(/3) = } 



j=0 



2- 



with 



1 



2sinh(2/3) 



(24) 



and c = {0, 0, 4, 16, 104, 416, 2 224, 8 896, 43 840, 
175 296, 825 648, 3 300480, 15101920, ...} up to order 323 
(at T = 1.5 the last term contributes « 0.28 x 10~ 158 ). 
The volume of the Wulff plot is given by [l^ | 



w 



where 



4 



dx cosh 



tr = 2 



cosh^(2/3) 



sinh(2/3) 



cosh(x) 



ln[tanh(/?)] 



(25) 



(26) 



is the interface tension of the (1,0) surface (i.e., in direc- 
tion of the axis ). F or the (1,1) surface the exact expres- 
sions reads (2lL 1221 



a x =^ln[sinh(2/3)] 
P 



(27) 
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B. Parameters for the n.n. Ising model on a 
triangular lattice 



Parameters for the n.n.n. Ising model on a 
square lattice 



The critical temperature of the triangular lattice is [2c 



T c = 



In 3 ' 



(28) 



For the spontaneous magnetisation Potts [2J] gave in 
1952 the expression 



mo = Wl- 



16exp(-12/3) 



[1 - exp(-4/3)][l + 3exp(-4/3)] 



(29) 



In contrast to the large number of low-temperature series 
expansions for the square lattice, we are aware of only 
two published papers for the triangular lattice [25], l26j |. 
In the second paper two more coefficients for the same 
series are given: 



xifi) = /3 2_\, Ciul u = ex P( — 4/3) 



(30) 



For the next-nearest neighbour model none of our pa- 
rameters are known exactly. The inverse critical tem- 
perature was given by Nightingale and Blote [28| using 
a transfer-matrix technique they call "phenomenological 
renormalisation" to be 



/3 C = 0.190192 69(5) . 



(35) 



In [29j this value was independently established using 
Monte Carlo simulations and finite-size scaling proce- 
dures. All other quantities are unknown in the literature 
and, therefore, computer simulations must provide the 
values. In the case of the magnetisation and the mag- 
netic susceptibility this is quite easy. A simple Monte 
Carlo algorithm at the desired temperature gives a time 
series of the magnetisation M. Then, the spontaneous 
magnetisation and the susceptibility are given by 



m Q 



1 N 

4^ 



(36) 



where c = { 0, 0, 4, 0, 48, 16, 516, 288, 
5 328, 3 840, 53 676, 45 488, 531600, 505 584, 5199 404, 
5 399 136, 50 369 760, 56 095 776, 484 296 732, 571 273 344, 
4628 107216 }. Finally, for the volume of the Wulff plot 
no explicit solution is available. Shneidman and Zia (27| 
showed the correct solution to be the integral 



and 



tt/6 

W(/3)=6 J d9 r 2 {6) 
o 

with a function r(9) given implicitly by 



(31) 



cosh 



3 + exp(2/3) 
-2 + 2exp(2/3) 

+ cosh [r/3 sin (9)] + cosh 



/3sin(| 



7T 

r/3 sin ( - 



(32) 



For the angles 9\ = I = 0,1,..., 11 the interface 

tension in direction normal to the equilibrium surface is 
given by r{0{). In the direction 9 = 7r/6 the minimal 
radius r m i n can be found to have the value 

fmin = "o = 



■ ( 33 ) 



The maximal radius r max is located at 9 — and Eq. 
simplifies greatly to 



2 (e^-\ 
r ma x = <J\ = — — In 



V3/3 



(34) 



= L 

X y 



i=l \ t=l / 



(37) 



where N is the number of Monte Carlo measurements 
and V = L x L the volume of the system. In the desired 
temperature range T ~ (2/3)T c the spatial correlation 
length ^ is very small and therefore already for moderate 
lattice sizes rather precise estimates can be achieved (4|| . 
Figure [6] shows the results of a Metropolis simulation of 
the n.n.n. Ising model at T = 4.0. 

To obtain the Wulff free energy is a much more de- 
manding task. Several methods are known, e.g. thermo- 
dynamic integration (30l . [3l| . Here, we will discuss two 
different ideas, namely a fit to the distribution of P(M) 
and a simple argument that the value of tw does not 
differ much from the appropriately scaled planar surface 
tension ctq. 

For our first method we exploit the fact that the prob- 
ability distribution for the largest droplet can be written 
as [32j 

P d cc exp (-/3rw\Ad) ■ (38) 

Using Eq. © and under the assumption Vd w vl the free 
energy in the exponent is 



(39) 



The assumption that the total overturned volume vl 
is consumed by the droplet volume v,i is certainly ful- 
filled the better the larger the droplet is. As is well 
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FIG. 6: The horizontal (green) line marks the mean of (a) 
the spontaneous magnetisation mo(L) and (b) the magnetic 
susceptibility x(L) f° r system sizes L = 40 . . . 1280 at T = 4.0 
of a n.n.n. Ising model. Its value gives an estimate for mo and 
X at L — > oo. Here, we read of the values mo = 0.947 2825(2) 
and x = 0.044 676(2). 
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FIG. 7: (a) Fit of the distribution In P d (M) = 
-/3rwv /l / 2 (l - M/M ) for a V = 160 x 160 n.n.n. Ising 
model at the temperature T = 4.0 in the range m — [0.4, 0.4 + 
400/160 2 ]. (b) Fit of the Wulff free energy tw vs. the inverse 
system size L at temperature T = 4.0 for L = 40, 80, . . . , 640. 
The error bars are obtained from (at least) 10 independent 
simulations per data point. 



known, the droplet can grow until it reaches the so-called 
droplet /strip transition point which is roughly located at 

M ds = M (l-^] . (40) 

With Eqs. fl35J) and (J3SJI, a linear fit of the form y = 
Twx + c can be achieved, where y — logPd and x — 
-/Vl/ 2 (! -M/Mq). Figure (a) shows such a fit for 
the 160 x 160 n.n.n. Ising model at the temperature 
T = 4.0 and for a range m = [0.4000,0.4156] which 
is close to the droplet/strip transition point located at 
TOds = Tno (1 — 2/V) w 0.3442. The data stems from 
a constrained multimagnetic simulation. To extract the 
value of the Wulff free energy in the thermodynamic limit 
of large systems, several simulations at different lattice 
sizes must be performed. In Fig. [7] (b) the scaling of the 
Wulff free energy is shown in dependence of the inverse 
lattice size. The intersection of the linear fit with the 
ordinate gives an estimate of ryy = 7.50 ± 0.02. 

Finally, we want to make three remarks about the 
given method. Firstly, we are fully aware of the fact, 
that Eqs. (|20|) and ([8|) give a "correction" to the fit done 
last. Using Vd(X) the fit would be valid for any droplet 



size up to the condensation/evaporation point and not 
only for large droplets nearby the droplet /strip transi- 
tion point. But on the other hand, the fit would not 
be a linear anymore and more important, the theoreti- 
cal predictions that we want to compare with would mix 
up with the parameter estimation. Secondly, it is pos- 
sible to measure during the simulation the droplet size 
i>d and fit directly TyvyA'd ms tead of P(m). Here, the 
disadvantage lies in the computational effort to measure 
the droplet size. While the magnetisation comes at no 
additional cost, a single measurement of the volume of 
the largest droplet needs 0(V) operations. Thirdly, we 
want to emphasise the importance of the initial starting 
conditions of the simulation. An ordered start where the 
first n spins point in one direction and the next V — n 
in the other direction is in fact a strip configuration. As 
discussed in 0, EB] between the strip configuration and 
the droplet configuration there is an exponentially large 
barrier that might not be overcome during the equilibra- 
tion phase, even so a droplet configuration has a much 
lower free energy for the constrained magnetisation range 
chosen. 

The second method to obtain r\v is based on the as- 
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sumption that, at the considered temperature, the in- 
terface tension for different angles 9 is roughly isotropic. 
This can be verified in detail for the n.n. Ising model, 
where the interface tension for an arbitrary angle 6 is 
known analytically [33|. For the planar interface the ex- 
pression (also given by Onsager [2^, 0, is a\ 



sq 

j 



2 J + Tin [tanh( J/T)\ and the expression for the "worst 
case", i.e. along the main diagonal of the lattice, is 
cr^ q = v^T In sinh(2 J/T) (also given by Fisher and Fer- 
dinand [H|). For all temperatures larger than T = 1.5, 
the relative difference of <7g q and cr^ q is smaller than 
1.3%. Obviously, the Wulff shape is still rather circu- 
lar at low temperatures and the quadratic form becomes 
only apparent close to T = 0. With this crude heuris- 
tics, the interface tension per unit volume at T = 1.5 

is 20Fc7o q = 4 - 219 - This is q uite close (99.37%) to the 
correct value = 4.245. An even better approximation 
is 2y/n(a s q + cr 1 q )/2 that deviates only 0.006% from the 
actual value. The same holds true for the triangular lat- 
tice. Using Eq. ||33J) one finds at T = 2.4 « §T C a value 
of 20Fcr* ri = 7.50657 which is only 0.005% smaller than 
the exact value of r^ 1 . Including Eq. (|34| for the im- 
proved estimation 2y / 7r(<7Q rl + a\ n )/2 yields a remarkably 
small difference of tiny 6 x 10~ 7 % to the exact result. A 
more detailed discussion concerning the approximation 
of a(0) can be found in [27| . For the n.n.n. Ising droplet 
the low-temperature Wulff shape is an octagon, i.e. it is 
much closer to the high temperature (low interface ten- 
sion) form, namely a circle. Therefore, it is reasonable 
to assume that above approximation might work as well. 
The planar interface tension can be measured using a 
multimagnetical (flat in the distribution of the magneti- 
sation) simulation, the result of which is a double-peaked 
magnetisation density P(m). In the limit of large system 
sizes L, it holds in two dimensions 36] 



In 



P 



(L) 



2(3a L 



(41) 



where P^l is the value of the density in the mixed phase 

region to rs and PmJ* the value at its maxima (m = 
imo). Figure[5] (a) shows the result of 13 multimagnetic 
simulations for the systems sizes L — 6 to L = 30 [29| . 
For every system the maximum and minimum probability 
-Pmax and P^l were read off and repeating the simula- 
tions ten times error bars were obtained. For L > 10 
the resulting values are plotted in Fig. [8] (b). An in- 
finite system size extrapolation in 1/L yields a value of 
do = 2.136±0.001 for the planar interface tension. Then, 
the estimate for the Wulff free energy (assuming a circu- 
lar droplet shape) is r w « x 2.136 = 7.571 ± 0.004 
which in fact is a lower bound, as the interface tension 
gets minimal along the directions of the interactions. 

Table [J gives the numerical values for the spontaneous 
magnetisation too, the susceptibility x and the Wulff free 
energy at the temperature T were the simulation took 
place. The temperature was chosen to be T « 0.66 T c - 
0.76 T c which is a good compromise between simulation 




(b) 



lit of simulation data — 

0.04 0.06 0.08 0.1 

inverse system size 1/L 



FIG. 8: (a) Distribution of the magnetisation m for the n.n.n. 
Ising model at T = 4.0 and system sizes L — 6,8, ... , 30. (b) 
Scaling of the interface-tension estimates from the histogram 
method: The straight line shows the fit bi^P^^c/P^J/L = 
2/?<7o(l + a /L) for L > 10 with goodness-of-fit parameter 
X 2 /d.o.f. = 1.1, yielding an planar interface tension estimate 
of ao = 2.136 ±0.001. 



speed (freezing at low temperatures) and compactness of 
the droplet (see the r.h.s. of Fig. [3] for a typical configu- 
ration) . 



D. Correction of the units in the parameter A 

After all constants are known, there are still some con- 
siderations to be made, before the parameter A can be 
calculated. The magnetisation too and the susceptibil- 
ity x are intensive quantities that follow from the cor- 
responding extensive quantities normalised (divided) by 
the volume. It is convention that for spin systems the 
volume is expressed by the number of spins, i.e. every 
spins accounts for a unit volume. In contrast, the free 
energy of the Wulff droplet is measured (again by con- 
vention) in units of the cell volume that is calculated 
given the lattice spacing a as input. As possible way to 
treat this situation is to normalise all quantities to cell 
volume, which would mean, that too and x are given in 
very unfamiliar units. We refrain from this step in or- 
der to keep things comparable to literature and instead 
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TABLE I: Numerical values for the magnetisation mo, suscep- 
tibility x an d Wulff interfacial free energy density rw entering 
the parameters A = A(vl, mo, Xj r w) defined in Eqs. (|23[) to 
(|36|) at the simulation temperature T for the three models 



studied. 




n.n. sq. 


n.n. tri. 


n.n.n. sq. 


T c 


2.269 


3.641 


5.258 


T 


1.500 


2.400 


4.000 


T/T c 


0.6610 


0.6592 


0.7608" 


mo 


0.9865 


0.9829 


0.9473 


X 


0.02708 


0.01959 


0.04467 


rw 


4.245 


7.507 


7.502 


2mo/rwX 


16.93 


13.14 


5.307 



"The temperature T = 4.0 was chosen without the knowledge of 
the critical temperature, certainly a value of T = 3.5 would have 
been more appropriate. 



modify Eq. (114|) in a very slight way. In order to do so, 
we define a scaling parameter Ant where all parameters 
are consistent with the conventions from literature 



3/2 



Ant = 2 



XTw L 2 



(42) 



Here, vl is the number of spins of the largest droplet 
including overturned spins, L 2 is the total number of 
spins of the system and mo, x are the magnetisation and 
susceptibility normalised to the total number of spins. 
The normalisation of the Wulff free energy rw does not 
change as it is given in terms of the unit volume in lit- 
erature. Secondly we define A uv where all quantities are 
given in terms of the unit volume which is the intended 
meaning by Biskup et ai, 



A„ v = 2- 



IT, 



2 ^3/2 



XT W V 



In this representation f2 is the volume of the largest 
droplet, V the volume of the total system, and /iq and 
X are the magnetisation and susceptibility normalised 
to the volume of the total system. If vq is the Voronoi 
volume of one spin [37| (the volume of the Wigner-Seitz 
cell of one spin) measured in units compatible with rw, 
then it holds 



O 
V 

Mo 

X 



vlv , 
L 2 v , 
M M 

V 



mo 



v L 2 
PV (< M 2 ) - <M> 2 ) 




— [(m ) (mo) ) = - 



(44) 
(45) 

(46) 
(47) 

(48) 
(49) 











a 

4 
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> 
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FIG. 9: The Wigner-Seitz cell of the (a) n.n.n. and (b) trian- 
gular lattice. It contains only one lattice site and all points 
within the cell are closer to this point than to any other lattice 
site. The red lines indicate the construction principle using 
the normals to the connection of a lattice to its neighbours. 



Now, a geometric "correction factor" a from Ant to A v 
can be defined as 



oAii 



lit 



(50) 



(43) Using Eqs. (f42]) - (f50]) a can be expressed as 



Auv 

Am 



'Itw V _ V "o y 



1 



X . 



ilL 2 v v 1 



(51) 



To conclude, using the parameters from literature as 
given in Table HI the abscissa must not be scaled with 
A but rather with A/^/v^ where vo is the Voronoi vol- 
ume of one cell. 

For the square lattice the Voronoi volume that a spin 
occupies is 1 x 1 which makes the correction factor trans- 
parent. The same holds for the n.n.n. lattice that has 
(by incident) the same geometry as the n.n. lattice, see 
Fig. [9] In the case of the triangular lattice the Voronoi 
cell is a hexagon. Figure[5](b) displays the situation. If h 
denotes the half of the lattice side length a, then it holds 
a = 2h. Every hexagon is made up of 6 small equilateral 
triangles of side length b (dotted line). The height of such 
a triangle is h which is given by h = b\/3/2. It follows, 
that b = a/y3. Now, the volume of a hexagon is given 
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by 

«o eX =^ 2 = ^(^) 2 = ^ 2 - (52) 

Finally, for a — 1, the geometric factor a for the trian- 
gular lattice is 



«tri = 




E. Droplet measurement 

As mentioned at the beginning of Sec. IIII1 one of our 
primary goals was the determination of the volume of the 
largest droplet. A possible advancement would be to set 
up a multimagnetic simulation and measure every sweep 
or so the droplet volume. While this is certainly possible, 
it is not advisable, as the determination of the multimag- 
netic weight factors W(m) 1/P(m) alone is a demand- 
ing task and in the following analysis there is no use for 
them. Instead, we arranged several simulations at fixed 
magnetisation m (micromagnetic). Inserting Eq. ([3]) in 
(fl"4"]) and solving for M gives the relation between the 
parameter A and the magnetisation M 

Solving Eq. (J53J) for A yields 



MM) = {Vmo _ Mf" , (55) 

2xt w V 

which shows that a fixed magnetisation results in a fixed 
value A(M). Therefore, we actually selected for every 
lattice 38 reasonable values A^ = {0.00, 0.10, . . . 16}, with 
an emphasis on the vicinity of A c . Using Eq. (|55|) a 
set of corresponding magnetisation values Mi, usually 
non-integer values, was obtained. A subsequent rounding 
to the next allowed value of the magnetisation (AM = 
±2) gave the final values for the simulation. To take the 
influence of the rounding into account, Eq. (|54|) was used, 
resulting in a second set A of slighty shifted (cx 1/W) 
values Ai that correspond to the rounded magnetisation. 

To enforce the constraint of constant magnetisation 
we use a Kawasaki update scheme where an up-spin is 
exchanged with a down-spin. Since the total number of 
up- and down-spins does not change, the magnetisation 
keeps its value as well. This type of non-local Monte 
Carlo moves can be accelerated using a table storing the 
spins sorted according to their direction. Here, one sweep 
accounts for V spin exchange attempts. 

After every sweep our simulation determines the vol- 
ume of the second-largest cluster which is (per defini- 
tion) the volume v,i of the droplet. This is done in two 
steps. First a Hoshen-Kopelman [38| algorithm performs 



a complete cluster decomposition. Thereby spins that 
are connected in the sense that they share a bond and 
have like orientation become a unique number. Figure flOl 
shows the situation for a spin-field and n.n. interaction. 
The largest (partially drawn) cluster (red) having clus- 
ter index 1 is the background, the cluster in the center 
(green) with cluster index 2 is the droplet we are look- 
ing for. Inside this droplet are smaller clusters located 
with cluster index 3, 4 and 5 (light blu e, y ellow, pur- 
ple). In the next step a flood-fill routine [39j . essentially 
a geometric depth first search, scans the droplet. Start- 
ing from an arbitrary position (that was recorded during 
the cluster identification step) it stops only when it finds 
spins that belong to the largest cluster (background). 
Thereby spins/clusters of opposite sign that lie within 
the droplet are subsumed. The result of this operation 
is shown in Fig. [10] (b). The thick blue line indicates the 
border between the droplet, i.e., cluster number 2 and 
all clusters which do not have the cluster number of the 
background, and the background. While this method is 
easy to implement and for the n.n. square lattice fool- 
proof, in case of the n.n.n. square lattice there are some 
pathological cases. Figure [11] shows such an ambiguous 
situation. FigurefTTI (a) presents the droplet as identified 
by our algorithm. In contrast, Fig. [11] (b) is an (imag- 
inary) alternative version resulting from the closing of 
the inclusion of background spins. The justification of 
the right pictures is given by the fact that the n.n.n. 
model has an interaction along the diagonal which con- 
nects the two surface spins (yellow). Fortunately, it is 
not necessary to decide upon which scenario is the more 
physical one. Every inclusion of reasonable size causes a 
large number of broken bonds due its surface. Therefore, 
configurations with inclusions are highly suppressed for 
temperatures well below the Curie point. To be on the 
safe side, we analysed several simulations of the n.n.n. 
square lattice for different system sizes with both meth- 
ods at the same time, i.e., for identical configurations 
the droplet was measured a second time with an algo- 
rithm that closes inclusions, to find negligibe differences. 
In the end we deciced to keep things as simple as pos- 
sible and therefore used only the combination Hoshen- 
Kopclman/flood-fill for our data generation. 



IV. NUMERICAL RESULTS 

For all three systems and at every value of A we per- 
formed simulations at five different lattices sizes L = 
40, 80, 160, 320, and 640. Every simulation ran at least 
20 000 sweeps for the thermalisation and at least 200 000 
sweeps for the measurements. To obtain the error bars 
reliably, 10 independent simulations were run for each 
data point. For the creation of pse udo random numbers 
we use the R250/521 generator [H El| . 

Having the numerical values of too, x, and r\y in place 
(see Sec. Mil) , the region of interest can be estimated. For 
A = 0.92 ss A c and the values from Table [I] correspond- 



10 





11 

T 






1 

T 




i 

T 


1 




1 






1 










t 


t 










1 

T 






t 


t 














{ 






t 


+ 


t 


t 










1 

T 






t 


41 














1 

T 






5 l 




t 


t 










1 










t 


T 










i 










T 


T 


















— 


1 


1 




























I 










1 


1 










i 

T 








li 


t 


t 










1 






t 


t 


1 












I 




t 


t 


t 


t 


t 










1 

T 




J 


J 


j 




t 










1 
T 










T 


* 

T 








_ 


T 
I 










t 


t 
















1 


1 








- 












I 


1 









FIG. 10: Cut out of a spin field, the red background cluster 
should be much larger, cf. the r.h.s. of Fig. [3] (a) The colors 
and the (small) numbers indicate the clusters detected and 
enumerated by the Hoshen-Kopelman routine, (b) The thick 
blue line surrounds the droplet (second largest cluster) found 
by the flood-fill routine. 
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FIG. 11: (a) Cut out of a spin field for the n.n.n. Ising model 
with a droplet (green) detected by the flood-fill algorithm. 
Apparently, the inclusion on the lower right side of the droplet 
(two spins) has a connection to the backgrund and does not 
count to the volume of the droplet, (b) Another way to in- 
terpret the situation where the spins are part of the droplet. 



ing to the n.n. Ising model, for L = 640 the magnetisa- 
tion is estimated with Eq. (|54|) to be m w 0.9827. To see 
the relevance of this figure we performed a multimagnetic 
simulation coupled with the parallel tempering algorithm 
[42| for the n.n. Ising model, the result of which can be 
seen Fig. [2j It shows the upper part (in the vicinity of the 
magnetisation peak in Fig. [T]) of the distribution of the 
magnetisation P(m) that exhibits for larger lattice sizes 
a clear cusp which divides the evaporated and condensed 
region. Within the evaporated region it has a Gaussian 
form according to Eq. l[7]). while in the condensed region 
a stretched exponential behavior is visible, cf. Eq. Q. To 
verify this quantitatively, Fig. [T^] shows a fit of a Gaus- 
sian curve and a stretched exponential curve to the upper 
part of the distribution of the magnetisation lnP(m) for 
the n.n. Ising model. The point of intersection m x is 
given by the condition 

Ve-^+ri = - (mx :" W)2 (56) 

the solution of which is a fourth order equation. With 
the parameters from the fit m max = 0.9864, a 1 = 1.042 x 
10~ 7 , c = 0.9858 and h — —1340, it evaluates to m x = 
0.9829, which is quite close to the aforementioned value 
calculated with Eq. (|54[) . The Gaussian fit which corre- 



sponds to the pure fluctuations part where A = can be 
compared to —j3Ff from Eq. 0. It yields for the suscep- 
tibility x = PVa 2 = 0.6666 x640 2 x 1.042 xl0- 7 « 0.028, 
a value quite close to the infinite-volume value given in 
Table H] of 0.02708. In the droplet dominated regime we 
have approximated the full mixed phase expression by 
neglecting the contributions of the fluctuations, which 
corresponds to putting A = 1 in Eq. ([TTj) . Over the fit 
range the neglected part contributes less than 4%. Even 
in the worst case, located at the cusp where A = 2/3, it 
amounts only to a value of approximately 9%. To ob- 
tain these values the ratio F d (l)/F(X) = 4n/A/(3A + 1) is 
evaluated using Eq. (f2"0"|) in conjunction with Eq. Q55p 
which yields an expression A = A(M). This is cor- 
roborated by the fact that, when fitting the droplet 
regime without fluctuations, from — j3F^(l) the Wulff 
free energy is approximated as t\v = —h/\j3y/V/ (2c)] = 
1340/[0.6666 x A /640 2 /(2 x 0.9858)] w 4.410, which is, 
again, quite close to the value of 4.245 given in Table fl] 

To have another "visual proof" that something differ- 
ent is happening on the two sides of the cusp in Fig. [2] we 
took several snapshots of the configurations that occurred 
during a simulation run. The two plots of Fig. [3] display 
an evaporated (left) and a condensed system (right), re- 
spectively. Both systems have the same number of over- 
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FIG. 12: Gaussian fit and stretched exponential fit of the 
distribution of the magnetisation P(m) for a the result of a 
L = 640 n.n. Ising simulation at T — 1.5. The left vertical line 
(magenta) indicates the transition magnetisation M(A C )/V 
predicted by Eq. (|54|1 while the the right vertical line (purple) 
coincides with the intersection point of the two fits. 



turned spins, i.e. the same magnetisation, which was cho- 
sen to be right at the transition point. While both config- 
urations occured during an actual simulation run, they 
do present extreme cases. When looked at the set of 
the largest cluster sizes recorded in the simulation run, 
the evaporated cluster configuration corresponds to the 
smallest number in the set and the condensed configura- 
tion corresponds to the largest number in the set. 

A final affirmation that the point under consideration 
was chosen correctly, can be derived from a look at the 
time series of the magnetisation m in Fig. [T2J The di- 
rect comparison shows a block structure in the time se- 
ries that coincides with cusp in the distribution P(m). 
Clearly, a sign for a barrier in the free energy. 

In Figs. Q3]- EI] w e show our main results, the frac- 
tion A(A) for the three observed lattices. The (black) 
solid line is the analytical value of A as shown in Fig. [5] 
Clearly, for larger lattice sizes the theoretical value is 
approached by the results of the simulation. Figure [TH 
(a) shows A in dependence of the magnetisation m. In 
Fig. Q3] (b) A is plotted for the same set of data points, 
but this time in dependence of A which essentially is a 
rescaling with v^ 2 . While in (a) the important region 
is barely visible, the rescaling leads to a blow up of the 
transition region making the theoretically predicted jump 
from Aa ~ to Aa ~ 2/3 at A c w 0.92 observable. This 
confirms that at the evaporation/condensation transition 
only 2/3 of the excess of the magnetisation goes into the 
droplet while the rest remains in the background fluctua- 
tions. The increase of Aa for A — * can be explained by 
the fact that the minimal cluster size is 1 and not an ar- 
bitrarily small fraction. In contrast, the excess that can 
be fixed analytically using Eq. (|14p can be much smaller 
than 1. 

In Fig. [T7] we compare A for L — 640 of the three 
different models. The nice agreement of the data points 
is a clear indication for the lattice independent universal 



behavior of the theory. An explanation for the slight 
discrepancy between the n.n. and the triangular lattice 
on the one side and the n.n.n. model on the other might 
be given by the slightly different temperature ratio T/T c 
(see Table HI. 



V. CONCLUSION 

Our Monte Carlo data clearly confirm the theoreti- 
cal considerations of Biskup et al. [f| Q for the case of 
the two-dimensional next-neighbour Ising system. While 
their results are only valid in the thermodynamic limit 
of large systems, we have shown that for practically ac- 
cessible sizes the theory can also applied. The observed 
finite-size scaling behavior fits perfectly with their pre- 
dictions for the infinite system. 

Moreover we have demonstrated that the theory, which 
to date has only been proven for the square lattice 
nearest-neighbour case, is actually universal in the sense 
that it is independent of the underlying lattice. The Ising 
model on the two-dimensional triangular lattice and on 
the two-dimensional next-nearest neighour lattice both 
approach the theoretically expected results nicely. Ap- 
parently, for the same relative temperature T/T c the 
finite-size behavior is identical. 

In order to achieve the correct scaling of the abscissa 
we presented several methods to estimate the Wulff free 
energy 7\y numerically. While in theory it should be 
straightforward to extract the value from the distribution 
of the magnetisation, due to limitations in the computer 
time for temperatures near the critical one, it can be more 
advantageous to resort to the isotropic approximation. 

All simulations were performed in thermal equilibrium 
and the abundance of droplets of intermediate size could 
be confirmed visually by looking at the distribution of 
droplets. We only state this fact here, while a more de- 
tailed analysis and the corresponding graphs will be pre- 
sented in a later publication together with more results 
on the finite-size scaling behavior of the systems and the 
shape of the free-energy barrier associated with the evap- 
oration/condensation transition. 
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FIG. 13: Time series of the magnetisation of a L = 160 n.n. Ising simulation at temperature T = 1.5. The distribution on the 
l.h.s. corresponds to the times series on the r.h.s and both were measured during the same simulation run. The lower (green) 
and upper (blue) horizontal lines indicate the transition magnetisation calculated from Eq. (I54[) for a lattice size L = 160 and 
L — > oo, respectively. The blocks in the time series are typical sign of a barrier in the free energy. 





FIG. 15: Fraction A for the two-dimensional triangular Ising 
model on square lattices of size L = 40, 80, . . . , 640 with pe- 
riodic boundary conditions at the temperature T = 2.4 « 

0.66 T c . Here, a = 1/Vua = Vg/V ^ « 1.075... is the ge- 
ometric factor, defined in Sec. IIII Dl The error bars are not 
plotted since their size is much smaller than that of the data 
symbols. The solid line shows the analytic solution in the 
limit L — > oo. 



FIG. 14: Fraction A for the two-dimensional n.n. Ising model 
on square lattices of size L = 40, 80, . . . , 640 with periodic 
boundary conditions at the temperature T = 1.5 « 0.66 T c . 
The error bars are not plotted since their size is much smaller 
than that of the data symbols. To show the influence of the 
scaling of the absissa, plot (a) and (b) use the same date. 
While in plot (a) the fraction A is given in units of the mag- 
netisation in plot (b) it is given in units of A. The solid line 
in plot (b) shows the analytic solution in the limit L — > oo. 
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FIG. 16: Fraction A for the two-dimensional n.n.n. Ising 
model on square lattices of size L = 40, 80, . . . , 640 with pe- 
riodic boundary conditions at the temperature T = 4.0 w 
0.76 T c . The error bars are not plotted since their size is 
much smaller than that of the data symbols. The solid line 
shows the analytic solution in the limit L — *• oo. 
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FIG. 17: Comparison of the fraction A for the three observerd 
Ising models (n.n., triangular and n.n.n.) for the size L = 640 
and the temperatures T = 1.5, 2.4, 4.0. 
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